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ABSTRACT 

We present the results of a numerical study on the effects of metal enrichment and 
metal cooling on galaxy formation and cosmic star formation (SF) history using cos- 
mological hydrodynamic simulations. We find following differences in the simulation 
with metal cooling when compared to the run without it: (1) the cosmic star formation 
rate (SFR) is enhanced by about 50 & 20% at z = 1 & 3, respectively; (2) the gas mass 
fraction in galaxies is lower; (3) the total baryonic mass function (gas + star) at z — 3 
does not differ significantly, but shows an increase in the number of relatively massive 
galaxies at z = 1; (4) the baryonic mass fraction of intergalactic medium (IGM) is 
reduced at z < 3 due to more efficient cooling and gas accretion onto galaxies. Our 
results suggest that the metal cooling enhances the galaxy growth by two different 
mechanisms: (1) increase of SF efficiency in the local interstellar medium (ISM), and 
(2) increase of IGM accretion onto galaxies. The former process is effective throughout 
most of the cosmic history, while the latter is effective only at z < 3 when the IGM is 
sufficiently enriched by metals owing to feedback. 

Key words: method : numerical — galaxies : evolution — galaxies : formation — 
galaxies : high redshift — galaxies : mass function — cosmology : theory 



1 INTRODUCTION 

Gas cooling plays a key role in galaxy formation. Accord- 
ing to the hierarchical structure formation model, galax- 
ies form in gravitationally collapsed dark matter halos. 
The dissipative baryonic matter cools by e mitting radiation 
and condenses into halos to form galaxies dRees fc Ostrikeri 
ll977l : IWhite fc Reedll978h . Subsequently stars form through 
the radiative cooling o f interstellar medium (ISM; e.g., 
iMcKee fc Ostrikerll2007l ). 

The cooling rate of gas varies with its chemical abun- 
dance, because different atoms and molecules have different 
cooling rates. The chemical abundance of cosmic gas evolves 
as a function of time and environment. In the early universe, 
the primordial gas consists of mainly hydrogen (mass frac- 
tion of X ~ 0.74) and helium (mass fraction of Y ~ 0.25) 
with a negligible amount of light metals (e.g. iBurles et al.l 
2001). Throughout this paper, all the elements heavier than 
helium are collectively called 'metals'. Most metals are pro- 
duced inside stars and spread out to the space by galactic 
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winds driven by supernova (SN) explosions. In the metal- 
enriched gas, the total cooling rate is significantly enhanced 
relative to that of the primordial gas by the atomic emission 
lines from recombination of ionised metals. 

Metals may influence galaxy formation in two differ- 
ent ways. Firstly, the SF efficiency is increased owing to 
the shorter gas cooling time. Star formation generally takes 
place in dusty metal-enriched ISM, therefore, the enhanced 
cooling rate by metals would boost up the SF efficiency. 
However, the energy /momentum feedback by SNe may sup- 
press the subsequent star formation after the initial star- 
burst by heating up the ambient gas. This complicates the 
situation, and it is not clear whether the net effect of metal 
cooling and SN feedback would be negative or positive. The 
feedback process is highly nonlinear, therefore, a direct nu- 
merical simulation would be a useful tool to explore the 
effects of feedback and metal cooling. 

Secondly, the metals dispersed into the IGM by SN feed- 
back also enhance the cooling of IGM. This can lead to the 
increase of IGM accretion onto galaxies, because colder gas 
can sink into galaxies more easily. Both of the above two 
effects may affect the galaxy growth and cosmic SF history 
significantly. Therefore it is essential to include the effects of 
metal cooling and chemical enrichment by feedback in the 
studies of galaxy formation and evolution, and capture the 
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complex two-way interactions between galaxies and IGM. In 
addition, metal enrichment also affects the equation of state 
by altering the mean molecular weight of gas, which influ- 
ences the hydrodynamic calculation of gas thermal state. 

Cosmological hydrodynamic simulations are widely 
used in the studies of gal a xy formation and cos- 



mic SF histo r y (e . g., ICen] fl99l IKatz et ail 1 19921: 



Cen fc Ostrikeil Il993l ; IKatz et all [ 



2000; iPearce et all [2001; 



1996 



Nagamine et al.l 
iNagamine et all l2001allbl; 



Ascasibar et al.l |2002|; iMurali et al.l |2002|; iNagaminel 
20021; IWeinberg et al.1 |2002|; iKawata fc Gibson! 120031; 



Kravtsovl 12003: Marri & White] 12003^ ISpringel fc Herriguistl 
2003bl: [Governato et all 120041; iRobertson et all 12004 
Weinberg et aJj 12004 INagamine et all 12004 l2005a~ 



Cen et alT 



Hocft et al 
20061; 



20051; iKeres et all 120051; IFinlator et al.1 12006 
120061: INagamine et al.ll2006l; IScannapieco et al.l 

Oppenheimerl 20071 



Tasker & Bryan 2006: Dave 



Finlator et al. 2007; Governato et al. 2007; Kobayashi et al. 
20071 ; iTornatore et al.ll2007l : IFinlator fc Davell2008l ') Most of 
these works included the treatments of radiative cooling by 
H & He, star formation, and SN feedback. However, some of 
the simulations did not include the effects of metal cooling 
and/or chemical enrichment by galactic wind. Furthermore, 
the effects of metal cooling on galaxy growth and cosmic SF 
history has not been explored in detail using cosmological 
hydrodynamic simulations and presented in the literature, 
as it is costly to run large cosmological simulations with 
and without the effect of metal cooling. 

In this paper, we investigate the effects of metal cool- 
ing and metal enrichment on galaxy growth and cosmic SF 
history using a series of cosmological hydrodynamic simula- 
tions with and without metal cooling. The aim of this paper 
is to single out the effects of metal cooling among our simula- 
tions. This work can be regarded as our initial step towards 
the long-term goal of developing more complete cosmologi- 
cal hydrodynamic code with physically motivated models of 
star formation and feedback. 

The paper is organised as follows. In § [51 we describe 
our simulation method focusing on how we implement the 
metal cooling. In § [3] we study the metal cooling effects on 
the global properties, such as cosmic SFR, phase space (p 
vs. T) distribution of gas, and the evolution of four phases 
(hot, warm-hot, diffuse, and condensed) of baryons. We then 
study the galaxy mass functions (§13]) and gas mass fractions 
(§0- Finally, we discuss and summarise our findings in §|6] 



2 NUMERICAL TECHNIQUE 

2.1 Simulation and Metal Cooling 

We use the updated version of the Tree-particle-mesh 
(TreePM) sm oothed particle hydrodynamics (SPH) code 
GADGET-2 (|Springell l2005h for our cosmological simu- 
lations. The gravitational dynamics is computed by a 
TreePM algorithm, which uses a particle-mesh method 
l|Hocknev fc Eastwood! Il988h for the long-range gravita- 
tional force a nd a Tree method fo r the short-range gravi- 
tational force ([Barnes fc Hut|[l986l ). This hybrid algorithm 
makes the gravitational force calculation faster than a 
Tree method and allows better force resolution than a PM 
method in dense regions. The gas dynamics is computed by 



an SPH method. The SPH is particularly useful if the sim- 
ulation needs to resolve large dynamical range, which is an 
inevitable requirement for the study of galaxy formation in a 
cosmological context. Therefore a TreePM-SPH simulation 
can provide a fast and high resolution calculation for both 
gravitational dynamics and hydrodynamics. 

The GA DGET-2 code adopts the en tropy-conservative 
formulation ijSpringel fc Hernquistl |2002| ) , which alleviates 
the overcooling problem that previous SPH codes suffered 
from. Our basic simulations include radiative cooling and 
heating p rocesses for hydrogen and hel i um u sing a method 
similar to Kat z. Weinberg fc Hernquistl l|l996h . An external 
UV background radiati on is treated as a spatia lly uniform 



photoionising radiation ijHaardt fc Madau| 19961). and mod i- 
fied to match the Lya forest observations l|Dave et al.lll999l ). 
Implementing star formation and feedback from first prin- 
ciples is not feasible in current cosmological simulations, 
because the spatial & mass-scales of molecular clouds are 
not resolved. Star formation and SN feedback are repre- 
sented by the subgrid multip hase ISM model developed by 
ISpringel fc Hernquistl (|2003 a) . In the multiphase scheme, a 
single gas particle represents both hot and cold gas. The 
stars are formed in the cold portion when the density exceeds 
a given threshold, p t h, which is derived self-consistently 
within the multiphase ISM model. The related parameters, 
such as the normalisation of gas consumption time and evap- 
oration efficiency of cold gas, are set to sat isfy the empirical 
Kennicutt-Schmidt law (|Kennicuttlll998d lbh. The SN feed- 
back returns some fraction of the cold gas to the hot phase, 
and increases the thermal energy of the hot gas. 

As an extension to the multiphase ISM model, the sim- 
ulation include s a phenomenological model for SN-driven 
galactic wind l|Springel fc Hernquistl 12003a). The galactic 
wind is particularly important for distributing the metals 
produced by SNe into the IGM. We use the strong kine- 
matic wind with a velocity of 484 km s -1 . It has been shown 
that this model produces favourable results for the lumi- 
nosity function of Lym an-break galaxies at the bright-end 
l|Nagamine et aTl 120041) and the Hi column density dis tri- 



bution function (jNagamine. Springel fc Hernquistl l2004h at 
z = 3, when c ompared to the runs w i thout the wind. 

However, iDave fc Oppenheimerl (|2007l ) pointed out the 
problems of this galactic wind model by comparing with 
the observations of C IV absorption lines in quasar spec- 
tra, and suggested that the momentum-driven wind model 
ijMurrav et al.l l2005) is a more viable model, which can carry 
more metals with lower wind velocities. In this paper we 
choose not to modify our galactic wind model in order to 
single out the effects of metal cooling , and to allow direct 
compa risons to the previous works (|Nagamine et alj|2004l . 
200510]). Here we focus on the effects of metal cooling on 
galaxy growth, while IDave fc Oppenheimerl |2007l ) focused 
on the C iv statistics of the IGM. 

The metallicity of gas particles are also tracked by the 
code, assuming a closed box model for each gas particle. The 
yield (y = AM mcta i/AM gas ) of 0.02 is assumed. In princi- 
ple, there could be a time delay between SN explosions and 
chemical enrichment of the ambient gas. Unfortunately, cur- 
rent cosmological simulations do not have sufficient resolu- 
tion to track the detailed mixing process of metals. In our 
simulations, we ignore this time delay and assume an instan- 
taneous mixing within each gas particle. This assumption is 
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Figure 1. The relationship between metallicity parameters in our implementation of metal cooling. The four panels show the following 
quantities as a function of metallicity Z/Zq, where Zq is the solar metallicity. [Fe/H], hydrogen mass fraction X, helium mass fraction 
Y, and oxygen abundance 12 + log(0/H). Panel (d) represents how a-particles contribute to the metal evolution. 



reasonable, because the overall time step of the simulations 
is longer than the mixing time-scale on small scales. 

Ideally we would like to track individual metals and 
compute the cooling rate of each element at each time- 
step of the simulation. This approach has been used to 
study the evolution of individual galax ies IIRecchi et al 



-ft- 

2001 ; iLia et al.l |2002| ; iKawata fc Gibson! 120031 ; IKobavashi 



2004 ) . For cosmological simulations, it requires a large mem- 



ory to track individual metals and th e computation be- 
comes exceedingly expensive (but see IScannapieco et all 
120051 : iMartfnez -Serrano et al.ll2008l . for such an attempts). 
Our simulations collectively track the total mass fraction 
of metals (Z = m mcta i/m gas ) as the measure of metallic- 



ity. We obtain the chemical abundance and cooling rates 
for a gi ven Z and temperature by interpolating the table 
given in lSutherland fc Dopital (|l993l . hereafter SD93) for the 
standard collisional ionisation equilibrium model, and cre- 
ate a lookup table. Metal cooling is implemented in the full 
range of metallicity in the table from SD93, from [Fe/H]=- 
3 to [Fe/H] = l. In addition, the metal cooling outside of 
this range is computed by extrapolating the table. We re- 
quire that the extrapolation does not go lower than the pri- 
mordial cooling rate. Th e primordial abundance pattern is 
used where [Fe/H] < -1 (|Wheeler et al.lll989l ; [Bessell et all 
Il99lf) and the solar abundance pattern is used where [Fe/H] 
> (| Anders fc Grevessd fl989h . For -1 < [Fe/H] < 0, the 
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abundance pattern is computed by interpolating between 
the primordial and solar abundance patterns. Note that, 
[Fe/H] = log(np e /njj) - log(np e ^ / n HQ )> where n is the 
number density. 

Figure [1] shows the relationships between the metal pa- 
rameters ([Fe/H], X, Y, and 12 + log(0/H)) and the total 
metallicity Z. Because the abundance pattern is fixed for 
a given Z, [Fe/H] is a monotonically increasing function of 
Z. The main difference between the solar and the primor- 
dial abundance pattern is the oxygen or a-particle enhance- 
ment. In the early phase of metal enrichment, Type II SNe 
are the main source of metals and they tend to generate 
more a-particles than SN la. Therefore, it is expected that 
a low metallicity gas to have more a-particles than the high 
metallicity one. As mentioned above, we use two different 
chemical abundance pattern (at low- and high-end of the 
Z values) to model this evolution in chemical composition. 
We compute the oxygen abundance based on the Table 4 of 
SD93. In our code, we incorporate the metal cooling effect 
by adding the additional cooling contribution from metals 
on top of the primordial cooling rate, which is com puted in 
a similar fashion as described in iKatz et al.1 (Il996t ). 

Metal enrichment changes the mean molecular weight 
fi, which is needed to compute the temperature and inter- 
nal energy of gas. In order to compute p, we need to know 
the electron number density n e , which varies with metallic- 
ity. We obtain the values of n e from SD93 in the same way 
as we obtained the metal cooling rates. Our electron num- 
ber density estimate ignores the electrons from photoionised 
metals by the UV background radiation, but we expect that 
the number of electrons from ionised metals will be small. 

Varying the mean molecular weight p also influences 
th e SF threshold density pth in the multiphase ISM model 
of lSpringel fc Hernquist (|2003al ) . as well as the cold and hot 
gas fractions of multiphase gas particles. In this paper, we 
consider two different treatments of p t h- In the 'constant p t h' 
mo del, we take the value o f pth as in the original formulation 
bv lSpringel fc Hernquist] (|2003aT ) and do not modulate p t h 
by the varying p. In the 'varying pth' model, we take into 
account of the modulation of p t h by the change in p. 



2.2 Simulation Setup 

We use three series of simulations with varying box size and 
resolution. The effects of metal enrichment and metal cool- 
ing are tested by comparing the runs listed in Table [1] The 
effects of galactic wind are also tested by turning on and off 
the wind. The adopted cosmological parameters of all simu- 
lations are (fi m , H A , Q b , a$ , h) = (0.3 , 0.7 , 0.04 , 0.85 , 0.7), 
which mostly concur with the results of the Wilkinson M i- 
crowave Anisotropy Probe (WMAP1) (|Spergel et al.ll2003h . 

Our naming conventions for the simulations are as fol- 
lows. The first part of the run name denotes the particle 
number and box size used for the simulation series. The run 
with the extension of 'mc' (for 'metal cooling') adopts the 
'constant pth' model for metal cooling, the run with the ex- 
tension 'mv' (for 'metal, varying') adopts the 'varying pth 
model, the run with the extension 'nw' has no galactic wind, 
and the run without any extension uses the cooling rates for 
the primordial chemical composition (i.e., only H and He). 

The 'N216L10' series is the fiducial simulation set in our 



study. We implement one simulation without galactic wind 
('N216L10nw' run). The 'N216L10' r un has the same phys- 
ical mo dels as the 'Q4' run used in ISpringel fc Hernquistl 
l|2003bl 'l and lNagamine et al.1 (|2004T ) with no metal cooling. 

To examine the resolution effect, we implement lower 
resolution simulations, the 'N144L10' series. Owing to the 
comparably small box size and missing long wavelength per- 
turbations, we evolve the N144L10 and N216L10 series only 
down to z = 2.75. To evolve the simulation toward lower red- 
shifts, we need a simulation with a larger box size. Therefore 
we implement the 'N288L34' series, which has a larger box 
size than other series, but a lower resolution. This series is 
evolved down to z = 1. The runs within the same series use 
the identical initial condition, which makes the comparison 
more robust and enables halo-by-halo comparison if needed. 



2.3 Galaxy finding method 

In cosmological simulations, dark matter, gas, and stars are 
represented by particles, and these particles are Monte Carlo 
representations of matter distribution. In order to study 
galaxy formation using a cosmological simulation, we need a 
working definition for simulated galaxies. In this paper the 
simulated galaxies are defined as isolated groups of star and 
gas particles, identified by a simplified variant of the SUB- 
FIND algorithm developed bv lSpringel et al1l|200l[ ). In more 
detail, the code first computes a smoothed baryonic density 
field to identify candidate galaxies with high density peaks. 
The full extent of these groups are found by adding gas and 
star particles to the groups in the order of declining density. 
If all -/V m iJ3 nearest neighbour particles have lower densities, 
this group of particles is considered as a new group. If there 
is a denser neighbour, the particle is attached to the group 
to which its nearest denser neighbour already belongs to. If 
two nearest neighbours belong to different groups and one of 
them has less than iV m i n particles, these the two groups are 
merged. If two nearest neighbours belong to different groups 
and both of them has more than -/V m j n particles, the particles 
are attached to the larger group, leaving the other group in- 
tact. In addition, the gas particles in groups should be denser 
than O.Olpth- In this paper we are not concerned with groups 
without star particles. The fraction of such groups are very 
small, therefore they do not affect our conclusions. 



3 GLOBAL PROPERTIES 

3.1 Cosmic star formation history 

Figure [2] shows the effect of metal cooling and galactic wind 
on the cosmic SF history for the N216L10 series. The com- 
parison of N216L10mc and N216L10 runs clearly shows that 
metal cooling enhances the cos mic SFR throughout all red- 
shifts by 20 - 30%. Previously. iHernouist fc Springell (|2003l . 
Fig. 13) argued that the SFR density is enhanced by metal 
cooling mostly at lower redshifts and hardly at high-z, how- 
ever, our simulations suggest otherwise. (See § [6] for more 
discussion on this point.) 



1 -^min is the minimum number of gas and star particles that 
constitute one isolated group. In this paper we set iV m i n = 32. 
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Run Box size ^_ N p ^ itidm^ m gas ^ z C nd Metal cooling^ WindJ^ 



N144L10 


10.00 


2 x 144 3 


2.42 x 10 7 


3.72 x 10 6 


2.78 


2.75 


No (primordial) 


Yes 


N144L10mc 


10.00 


2 X 144 3 


2.42 X 10 7 


3.72 x 10 6 


2.78 


2.75 


Yes (const, pth) 


Yes 


N216L10 


10.00 


2 X 216 3 


7.16 x 10 6 


1.10 x 10 6 


1.85 


2.75 


No (primordial) 


Yes 


N216L10nw 


10.00 


2 x 216 3 


7.16 x 10 6 


1.10 x 10 6 


1.85 


2.75 


No (primordial) 


No 


N216L10mc 


10.00 


2 x 216 3 


7.16 x 10 6 


1.10 x 10 6 


1.85 


2.75 


Yes (const. p t h) 


Yes 


N216L10mv 


10.00 


2 x 216 3 


7.16 x 10 6 


1.10 x 10 6 


1.85 


2.75 


Yes (varying p th ) 


Yes 


N288L34 


33.75 


2 x 288 3 


1.16 x 10 s 


1.79 x 10 7 


4.69 


1.00 


No (primordial) 


Yes 


N288L34mc 


33.75 


2 x 288 3 


1.16 x 10 8 


1.79 x 10 7 


4.69 


1.00 


Yes (const. p th ) 


Yes 



Table 1. The simulations employed in this paper 



a Box size in units of h 1 Mpc. 

^ Initial particle number of dark matter and gas particles. 
c Dark matter particle mass in units of fi"'M0. 

^ Initial gas particle mass in units of /i -1 Mq. The star particle mass is set to 0.5m gas in our simulations. 

e Comoving gravitational softening length in units of /i _1 kpc. This is a proxy for the spatial resolution of the simulation, and the 
N216L10 series has a spatial resolution of physical 0.5 h~ 1 kpc at z = 3. 

J 'Primordial' means the cooling rate with primordial chemical composition. See ^ 12.11 for the description of 'constant p t h' and 'varying 
Pth' models. 

9 If employed, the wind is a 'strong' wind with v w = 484kms — 1 . 




Figure 2. Cosmic star formation histories for the N216L10 se- 
ries, comparing the runs with and without metal cooling. It also 
shows the run without a galactic wind (N216L10nw), and a lower 
resolution run (N144L10mc). Overall the SFR density is increased 
by ~ 20 — 30% by metal cooling throughout the entire redshift 
range. 



Figure 3. The same as Figure [2] but for the N288L34 series, 
which goes down to z = 1. The N216L10mc run is also shown for 
a resolution comparis on. The green dot-dashed lin e (HS03) is for 
the analytic model o f lHernauist fc Springell j2003l. Eq. (2)). 



Metallicity needs to be higher than ~ 10~ 2 Zq to no- 
ticeably increase the gas cooling rate, because the metal 
cooling rate is much smaller than that by H and He at 
Z < W~ 2 Z Q . The metallicity of diffuse IGM is generally 
lower than 10~ 2 Zq at z > 3, and the metal cooling effect 
on IGM is small. However the local SFR could be enhanced 
significantly even at high-z when the ISM is enriched to 
Z > W~ 2 Z e , which can be easily achieved after a few events 
of SN II. Our results suggest that star formation at high-z 



can also be enhanced by metal cooling, and not just at low-z. 
We also find that the peak of the cosmic SF history hardly 
shifts by the int roduction of metal coo li ng, co nsistently with 
the estimate bv lHernauist fc Springe] (|2003l ). 

The N216L10mv run has a lower SFR at z > 9 and 
a higher SFR at z < 7 than the N216L10mc run. This is 
because at high-z, the value of p t h is higher in the W 
run than in the 'mc' run due to lower metallicity and mean 
molecular weight, leading to a larger fraction of gas not being 



© 0000 RAS, MNRAS 000, 000-000 



6 Choi & Nagamine 



able to satisfy the SF criteria in the N216L10mv run than 
in the N216L10mc. The opposite is true at z < 3. 

Figure also shows that the SFR in the N216L10nw 
run is higher than other runs by an order of magnitude. 
This is because the gas can cool very efficiently without 
being ejected by the galactic wind. Although the current 
galactic wind models in cosmological simulations are not well 
established yet and need to be improved, both theoretical 
and observational studies evi dently suggest the important 
role of galactic wind feedback |Springel fc Hernquisdl2003al ; 
lMartinll2005l , l200rj ; lOppenheimer fc Davell2006l ). 

As expected, the N144L10mc run has lower SFR den- 
sities at z > 7 compared to the N216L10 series. This is 
because a lower resolution run is unable to resolve low-mass 
halos that would otherwise forms stars in them. When the 
star formation at high-z is underestimated, some of the gas 
remain unconverted into stars until lower redshift, resulting 
in a higher SFR at z < 4 in the N144L10mc run than in the 
N216L10mc run. 

Figure shows the cosmic SF history down to z = 1 
for the N288L34 series. Again, the SFR density in the 
N288L34mc run is systematically higher than that in the 
N288L34 run at all redshifts. Moreover, the difference in 
SFR density between the two runs increases at z < 3. By 
z ~ 3 the IGM metallicity becomes high enough to notice- 
ably increase the cooling rate, and the accretion of IGM onto 
galaxies becomes more efficient, which later fuels the star 
formation. Our results suggest that the cosmic SFR density 
is enhanced by metal cooling at z < 3 by both of the follow- 
ing two effects: 1) more efficient IGM accretion onto galaxies 
due to higher metallicity of IGM, and 2) increased SF effi- 
ciency due to higher cooling rate with the contribution from 
metal cooling. 

In Figure [3] the analytic model of Irlernquist fc Spri ngel 
l|2003l . Eq. (2), hereafter H&S model) is also shown in the 
green dot-dashed line. This model is a result of many cos- 
mological SPH simulations with increasing resolution, and 
the authors gave a theoretical interpretation to the empirical 
fitting formula for the cosmic SFR density. Since our runs 
adopt a lower value of <rg = 0.85 than the original value 
used to derive the H&S model (<7 8 = 0.90), we modify the 
H&S model by rescaling the j3 parameter in their Eq. (2) 
and Eq. (47). This reduces the SFR density at z > 5 by 
~ 10%, but not much at z < 5. The results of the N288L34 
and N216L10mc runs agree well with the H&S model at 
z < 4 and z < 6, respectively. This is expected, because 
the N288L34 run adopts the same physical models as in the 
original simulations that were used to derive the H&S model. 
At high-z, both runs fall short compared to the H&S model 
owing to insufficient resolution. 



3.2 Phase space distribution of cosmic gas 

The metal enrichment and metal cooling also change the 
phase space distribution (p — T diagram) of cosmic gas. Fig- 
ure|4]shows the state of cosmic gas at z = 3 for the N216L10 
series. As shown in panel (a), the baryons in the Universe 
can be broadly categorised into four different phases accord- 
ing to their overdensity a nd temperature: 'hot', 'warm-hot', 
'diffu s e', and 'condensed ' (|Dave et al.lll999l : ICen fc Ostrikerl 
1 19991 : iDave et all l200lf) . The first two phases are mostly 
the shock-heated gas in clusters and groups of galaxies. 



The 'diffuse' phase is mostly the photoionised IGM with 
lower temperature, which can be observed as the Lya for- 
est. The tight power-law relationship between p and T at 
P/P ~ 3 and T < 10 4 K is go verned by the ionisation equi- 
librium dHui fc Gnedinlll997h . where p is the mean density 
of baryons at z — 0. 

The 'condensed' phase is the high density, cold gas in 
galaxies. Because the primordial cooling curve has a sharp 
cutoff at T ~ 10 4 K, the temperature of condensed phase de- 
creases only down to this temperature in Figure [4ji,b. The 
characteristic finger-like feature at p/p > 10 4 is due to the 
pressurisation of the star- forming gas by SN feedback in the 
multiphase ISM model of ISpringel fc Hernquistl (|2003ah . A 
fitting f ormula for this effectiv e equation of state was de- 
rived bv lRobertson et al.l (|2004T ). which is shown by the blue 
dashed lines. Above the SF threshold density p t h, the gas 
particles in the simulation go into the multiphase mode and 
are allowed to form stars. 

The most noticeable and interesting change in the p — T 
plot by the metal cooling is in the distribution of 'condensed' 
gas. Star formation and SN explosions occur in the mul- 
tiphase star-forming gas, therefore the immediate effect of 
metal cooling appears in the condensed phase. The presence 
of gas at T < 10 4 K with metal cooling results from the ef- 
fective temperature calculation for the multiphase medium, 
and not from the direct gas cooling d own to T < 10 4 K. Ac- 
cording to the multiphase scheme in ISpringel fc Hernquistl 
(|2003al ). the star forming gas consists of two phases: hot gas 
(10 5 K < T < 10 s K) and cold gas (T = 10 3 K). The effective 
temperature for the multiphase gas is computed from the in- 
ternal energy of hot and cold gases, which are derived using 
the cold gas fraction. The cold gas fraction is computed from 
the cooling rate of hot gas and the gas density. The original 
GADGET imposes the effective temperature of 10 4 K at the 
SF threshold density. In our metal cooling scheme, however, 
the cooling rate of hot gas increases, and the cold gas frac- 
tion increases. Therefore, the multiphase gas above the SF 
threshold density with higher metallicity could have lower 
effective temperature of T < 10 4 K in our implementation. 

In the N216L10mc run, the value of pth was kept fixed 
as the original value, which causes a sharp edge at p/p — 10 4 
for the multiphase gas. In the N216L10mv run, the value of 
pth was varied according to the change in mean molecular 
weight, therefore the condensed gas has a spread in both 
temperature and density near p/p — 10 4 . 

Galactic wind returns the metal-enriched, star-forming 
gas to the hot and low density IGM, as evidenced in the p— Z 
diagrams. This metal-enriched, low-density (p/p < 10) gas 
is absent in the N216L10nw run, which clearly shows the ne- 
cessity of galactic wind to explain the metallicity observed in 
Lya forest. The phase distribution of IGM in the N216L10, 
N216L10mc, and N216L10mv runs are very similar. This is 
presumably because the effects of metal cooling is not so 
significant yet in the IGM at z > 3. As we will discuss in 
the next section, this is not be the case at z < 3. 



3.3 Evolution of four phases of baryons 

The connection between the four phases of cosmic gas 
is important in galaxy formation, and it is beneficial for 
us to study their evolution for a better understanding of 
metal cooling effects. A conventional view is that the dif- 
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(c) N216L10mc (d) N216L10mv 

Figure 4. The phase space and metallicity distribution of the cosmic gas at z = 3 for the four simulations of N216L10 series. For each 
run, three panels are shown for p — T, Z — T, and p — Z relationship. For plotting purpose, we plot only randomly selected 1% of the total 
gas particles in the simulatio n. The blue dashed line is based on the fitting function to the effective equation of state for the multiphase 
star-forming gas provided by Robertson ct al ] (|2004h . 



fuse gas is shock-heated to hot or warm- hot phase, which 
later cools down to t he condensed pha s e llDave et al.|[l999l ; 
ICen fc Ostriker1ll999l ; IPavI et al. 1 120011 ) . iKeres et al l (|2005l ) 
also suggested that some diffuse gas can migrate to the con- 
densed phase without being shock-heated, which they called 
the 'cold mode' accretion. 

Figure [5] shows the evolution of mass fractions of dif- 
ferent phases. All of our simulations show that the diffuse 



phase continuously decreases from high-z to low-z, while the 
'warm-hot + hot' and 'condensed' phases continuously in- 
crease, consistently with the general expectations. The IGM 
(i.e., hot + warm-hot + diffuse) is accreted onto galaxies, to 
become the condensed phase. The diffuse component (panel 
a) decreases from ~ 95% at z — 6 to w 55% at z = 1. 
The 'warm-hot + hot' component (panel b) increases from 
almost zero at z — 6 to ~ 35% at 2 = 1. The mass fraction 
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Figure 5. The evolution of mass fractions for the four phases of baryons: hot, warm-hot, diffuse, and condensed. Panel (b) shows the 
addition of the two phases, 'warm-hot + hot', and panel (c) shows the total masses in the IGM, i.e., 'diffuse + warm-hot + hot'. Here 
we compare the N216L10 and N288L34 series. 



of the total IGM (panel c) decreases from 96 — 98% at z = 6 
to 92 — 93% at z — 3. This reflects the increase of the con- 
densed phase (panel d) from a few percent at z = 6 to ~ 7% 
at 2 = 3. The se mass frac t ions a re in general agreement with 
the results of iDave et all (|200ll ). 

The comparison between N288L34 and N288L34mc 
runs shows interesting differences. At z > 3, the mass frac- 
tions of all four phases are similar in the two runs. However, 
at z < 3, the IGM mass fraction continues to decrease in the 
N288L10mc run, while it becomes almost constant in the 
N288L10 run. A similar behaviour is seen for the condensed 
phase, and the N288L10mc continues to increase its con- 
densed mass. This suggests that the metal cooling enhances 
the IGM accretion onto galaxies mostly at z < 3, which sup- 
plies the fuel for star formation. The cooling of warm+hot 
phase into diffuse phase is also enhanced, leading to a larger 
mass fraction of diffuse component in the N288L34mc run 
than in the N288L34 run. (The same is true for N216L10mc 
and N216L10 runs.) 

The N216L10nw run shows stark differences from all 
other runs. It has a significantly lower IGM mass fraction, 
and a very large mass fraction in the condensed phase. This 
shows that the IGM is mostly enriched by galactic wind, and 
that most of the masses are returned into the 'warm-hot + 
hot' phase, but not to the diffuse phase. 



4 GALAXY STELLAR & BARYONIC MASS 
FUNCTION 

4.1 At z > 3: 

Galaxy stellar mass function (GSMF) tells us how the stel- 
lar mass is distributed in different galaxies, and its redshift 
evolution reflects the growth of structure in the hierarchical 
universe. Figure [6] compares the GSMF in different simu- 
lations. We only show the range of M star > 2.5 x 10 7 M o , 
which corresponds to the limiting mass of 32 star particles in 
N216L10 runs. We chose 32 star particles as our resolution 
limits because all mass function start to turn over around 
this mass range. Panel (a) shows that the run with no galac- 
tic wind (N216L10nw) completely overestimates the GSMF, 
and the run with lower resolution (N144L10) underestimates 
the number of lower mass galaxies at M s t ar < 10 8 Mq rela- 
tive to the N216L10 run. Panel (6) shows that the metal cool- 
ing increases the masses of galaxies by ~ 20% at z — 3. As 
a result, both the N216L10mc and N216L10mv runs show a 
slight increase in the number of galaxies at the massive end, 
with a slightly stronger enhancement in the N216L10mv run. 

Figures |6j;,d demonstrate that the metal cooling also 
changes the redshift evolution of GSMF. In the N216L10mc 
run, the peak of mass function is shifting more toward the 
massive side than in the N216L10 run at z = 3. From z = 4 
to z = 3, the low-mass galaxies merge to form more massive 
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Figure 6. Stellar mass functions of simulated galaxies. Panel (a): For the runs without metal cooling at z = 3. The vertical dotted line 
indicates the galaxy masses with 32 star particles for N144L10 run. Panel (b) compares the N216L10 series at z = 3. The bottom two 
panels show the redshift evolution of GSMF in the N216L10 {panel c) and N216L10mc (panel d) run. 



galaxies, increasing the number of galaxies with M sta r > 
10 s M furthermore. 

It also seems that the formation of low-mass galaxies 
(Af s tar < 10 8 Mq) is enhanced by metal cooling at early 
times (z = 5 — 6). However we find that the behaviour of MF 
at the low-mass end is somewhat dependent on the threshold 
density for the grouping. As mentioned in § 12.31 our group- 
ing code imposes a threshold gas density of O.Olpth for a gas 
particle to be part of a galaxy. When we lower this thresh- 
old density, we find that the two MFs for the N216L10 and 
N216L10mc converged, as more particles are incorporated 
into galaxies at the outskirts of galaxies. Metal cooling in- 
creases the gas density in galaxies, enabling more low-mass 



galaxies to satisfy the gas density threshold and to be iden- 
tified as simulated galaxies. Therefore the behaviour of MF 
at M s tar % 5 x 10 7 Mq is somewhat dependent on the choice 
of the outer threshold density of galaxies. 

Figure [7] compares the gas and baryonic (star + gas) 
MFs for the N216L10 series at z — 3 (panels a & b), and the 
redshift evolution of the baryonic MF from z = 5 to z = 3 in 
the N216L10 (panel c) and N216L10mc (panel d) runs. Here 
we show only the range of M s tar > 5.0 x 10 7 M Q , which cor- 
responds to the limiting mass of 32 gas particles. Panel (a) 
shows that more gas is converted into stars in the runs with 
metal cooling (N216L10mc and N216L10mv) compared to 
the N216L10 run, while panel (b) shows that the baryonic 
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Figure 7. The stellar (panel a), gas (panel 6), and the total baryonic (panel c) mass functions for the N216L10 series at z = 3. 



MF is similar in all the runs. As we discussed in § 13.31 the 
IGM accretion onto galaxies is not so much enhanced yet 
before z ~ 3 due to relatively low IGM metallicity, which 
explains the similarity in the total baryonic MF. In the runs 
with metal cooling, the peak of GSMF is shifted toward 
higher mass, whereas the peak of gas MF is shifted toward 
lower mass. This suggests that the enhanced GSMF at z — 3 
is due to the increased SF efficiency by metal cooling within 
the galaxies. Again, the apparent enhancement in the num- 
ber of low-mass galaxies in the N216L10mc run compared to 
the N216L10 run is somewhat dependent on the threshold 
density of grouping, therefore it should be interpreted with 
caution. 



4.2 At 2 < 3: 

As the accretion of IGM onto galaxies become more efficient 
at z < 3 due to the increased IGM cooling by the metals, 
we expect that the total baryonic mass of galaxies would 
be more enhanced at z = 1 than at z = 3. Figure [8] shows 
the MFs at z = 3 and z = 1 in the N288L34 series, and it 
clearly demonstrates that this expectation is true. Here we 
show only the range of Af gas > 8.2 x 10 8 Mq, which corre- 
sponds to the limiting mass of 32 gas particles. The general 
trend in the three (star, gas, and total baryon) MFs is sim- 
ilar to that we saw in Figure [7] although the difference at 
2 = 3 between the N288L34 and N288L34mc runs is slightly 
smaller than in the N2I6LI0 series due to poorer resolution. 
Figure [8}' shows most prominently the enhancement of the 
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M [M ] 



Figure 8. The stellar, gas, and baryonic (star + gas) mass functions at z = 3 (left) and 2 = 1 (right) for the N288L34 series. 



total baryonic MF at z = 1 in the N288L34mc run owing 
to the metal cooling. In panels (e) & (/), the N288L34 run 
actually has a longer tail at the most massive-end than the 
N288L34mc run. The reason for this feature is not fully clear, 
but it may be related to the balance between IGM accretion 
and feedback. Owing to metal cooling, IGM accretion rate 
increases in the 'mc' run, and the SFR is also enhanced, 
leading to a stronger feedback. The amount of mass loss is 
greater in low mass galaxies, but the net heating of IGM 
is more significant for massive galaxies. The 'mc' run has 
stronger feedback, therefore its feedback heating may be- 
come more significant than IGM accretion for very massive 
galaxies. The significant IGM heating suppresses the growth 
of massive-end of mass function from z = 3 to z = 1 and 



results in shorter tail for the 'mc' run at z — 1 as shown in 
Figure ISM. 



5 GAS FRACTION 

For a fixed galaxy baryonic mass, the enhancement of star 
formation by metal cooling would decrease the gas mass frac- 
tion, /ga S = Mgas/Mbaryon- Figure [9] shows /ga S as a func- 
tion of galaxy stellar mass for the N216L10 series. Panel 
(a) shows the data only from N216L10mc run at z = 3, 
and each data point corresponds to a simulated galaxy. 
To characterise the distribution, we compute the following 
two quantities in each logarithmic stellar mass bin: 'aver- 
age' and 'median'. The 'average' is the ratio of total gas 
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Figure 9. The gas fraction of galaxies as a function of stellar mass. Each data point in panel (a) is a simulated galaxy in the N216L10mc 
run at z = 3. Here we overplot two lines representing the distribution: 'average' and 'median' (see text for details). Panels (ft) — (d) 
compare the gas fractions ('average' case) in the N216L10 series at z = 3, 4, Sz 5. 



mass to total baryonic mass for all the galaxies in each 
mass bin, i.e., J2i ^gas,i/ J2i Mbaryon,i- The 'median' case 
is simply the median of / gas values in each mass bin. Both 
quantities show a similar trend, however, there is a sharper 
drop-off at M star — 2.5 x lO 7 Af0 for the 'median' case 
in Figure [9^. This mass-scale corresponds to 32 star par- 
ticles in the N216L10 series. We find that there are many 
galaxies with / gas = above this mass-scale, which causes 
the sharp drop-off in the 'median' line. Below this limit- 
ing mass, galaxies are not resolved well, which results in 
an underestimate of star formation and an overestimate of 
/ gas . If we had a higher resolution simulation with finer 
particle masses, this limiting mass-scale would shift to a 



lower mass. Therefore the location of this sharp drop-off 
is currently determined by the resolution of our simula- 
tion. However, dark matter halos would stop forming stars 
at some lower limiting halo mass, if we had an infinitely 
high-resolution simulation, This lower limit to the galaxy 
mass is presumably determined by t he photoevaporation of 
gas by the UV background radiation | Recs 1986; Efst athioul 



1992j;iQuinn et alJl99d:lGnedirj200C 



Nagamine et al.ll2004l ; 



iPontzen et al. 2008; O kamoto et al 20081 ). Recent works 
suggest that star formation could be suppressed by the UV 
background in halos with Mhaio < 10 9 M@ at z ~ 3. Galax- 
ies with M s tar — 2 x 10 7 Mq would reside in halos with 
Mhaio ~ 2 x 10 9 Mm. The N216L10 series would resolve such 
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a halo with ~ 280 dark matter particles, and its mass res- 
olution is actually close to the astrophysical limit for dwarf 
galaxy formation at z ~ 3. Therefore the sharp drop-off in 
/ gas at Mstar — 2 x If) 7 Mq may not be so far from the true 
answer. 

We find that / gas increases with decreasing M sta r at all 
redshifts, regardless of metal cooling and wind effects. This 
trend is qualitat ively consistent with current observations. 
lErb et all (|20Q6Q estimate the gas fraction as a function of 
stellar mass using the rest-frame UV-selected star-forming 
galaxies at z ~ 2, and show that the gas fraction decreases 
with increasing stellar mass. And using the mean gas and 
stellar mass, they find the average / gas ~ 0.35, which agrees 
with the predicted gas fraction for massi ve galaxies in our 
simulation. In addition. iGeha et al.l (|2006T ) reported that the 
average neutral gas fraction is {/ ga s) = 0.6 for the local dwarf 
galaxies selected from the Sloan Digital Sky Survey. In our 
N216L10 series with metal cooling, the 'average' / gas reaches 
0.6 for galaxies with M star ~ 2 x 10 7 M at z = 3 & 4. 

Figures[9p,c,d show the redshift evolution of / gaa for the 
N216L10 series. In the runs with metal cooling (N216L10mc 
and N216L10mv), / gas is lower than in the N216L10 run 
by 20 — 30% at all redshifts for galaxies with M star = 
10 7,5 — 10 9 Mq. This result suggests that the metal cool- 
ing reduces / gaa owing to more efficient star formation. The 
values of / gas seem to be more convergent at the massive- 
end (Mat ar > 10 9 Mq). In addition, we find that galaxies in 
the N216L10nw run are the most gas-rich, because almost 
no galactic gas is returned to the IGM. 

Figures |9j3,c,d also show that / gas is higher at higher 
redshifts. For example, the well-resolved galaxies with 
Mstar = 10 9 M e have / gas ~ 0.45 at z = 5, but / gas = 0.25 
at z = 3. As we saw in Figure [BJi, there is a considerable 
growth in GSMF from z — 5 to z — 3. Together with the 
decrease in the gas fraction, these results suggest that most 
of the gas accreted during z = 3 — 5 has been converted into 
stars. 

Figure [TUb highlights the redshift evolution of / gas from 
z = 5 to z = 1 in the N288L34mc run. The gas mass fraction 
clearly decreases with decreasing redshift as the gas is con- 
verted into stars. The rate of decrease is greater at z — 5 — > 3 
than at z = 3 -> 1, with / gas ~ 0.5, 0.25, & 0.2 for z = 5, 3, 
& 1, respectively, for galaxies with M star > 3 x 10 9 Mq. 

Figure 110b shows the effect of resolution on / gaa . The 
location of the drop-off in / gas shifts to lower masses as the 
resolution is increased from N144L10mc to N216L10mc run. 
The location of the drop-off is at around our resolution limit 
- the vertical lines in Panel (b). Therefore the values of / gas 
are overestimated in unresolved galaxies. 



6 DISCUSSIONS AND CONCLUSIONS 

Using cosmological hydrodynamic simulations with metal 
enrichment and metal cooling, we studied their effects on 
galaxy growth and cosmic SFR. Owing to metal cooling, the 
SFR density increases about 20% at z = 3 and about 50% 
at z = 1. Our results suggest that metal cooling enhances 
the star formation through two different processes: 1) more 
efficient conversion of local ISM into stars (i.e., increase of 
local SF efficiency), and 2) the increase of IGM accretion 
onto galaxies. The former process is in effect essentially at 



all times at z < 15, because the local ISM can be instan- 
taneously enriched by SN explosions. This process enhances 
the SFR, but does not noticeably change the total baryonic 
mass of galaxies. The latter process, on the other hand, can 
increase the total baryonic mass of galaxies, as well as en- 
hancing the overall SFR density by supplying more gas for 
star formation. This process becomes effective only at lower 
redshifts (z < 3), because it takes some time to enrich the 
IGM up to Z > 10' 2 Z e . 

In this paper, we used two different prescriptions for 
star formation in the runs with metal cooling: 1) the 'con- 
stant /9 t h' scheme and 2) the 'varying pth' scheme. The value 
of p t h was fixed in the former scheme, while it was modu- 
lated according to the metallicity in the latter scheme. Both 
schemes show similar increase of cosmic SFR density (see 
Figure [SJ) , which implies that the overall SF enhancement 
by metal cooling is not so sensitive to the choice of p t h (but 
with some differences as we discussed in § I3.1J1 . 

The multiphase ISM m odel for star formation by 
ISpringel fc HernauistJ 12003a - ) contains two free parameters: 
Pth and the normalisation of gas consumption time-scale, t°. 
The values of these two parameters were origin ally chosen 
to ma tch the empirical Kennicutt-Schmidt law ( Ke nnicuttl 
Il998al lbh using simulations of isolated disk galaxies. Since 
we did not change the value of t+ in the runs with metal 
cooling (as well as p t h in the 'mc' runs), it is possible that 
our simulations may now violate the Kennicutt law. 

However, as we show in Figure [TT1 the plots of iVm vs. 
Ssfr for the N216L10 and N216L10mc runs are not so dif- 
ferent from each other. This can be understood as follows. 
As Figure[5j; showed, metal cooling increases the density and 
lowers the temperature of star-forming gas. The fundamen- 
tal scaling relationship between SFR and cold gas density 
does not need to change when the metal cooling is intro- 
duced; the star-forming gas particles would simply slide up- 
ward along the Kennicutt law, Esfr oc iVm 1 ' 4 . 

In the case of N216L10mv run, we varied p t h accord- 
ing to the gas metallicity, primarily scattering it to lower 
values because higher metallicity increases the cooling rate 
and lowers p t h- Therefore more gas particles become eligible 
to form stars and Esfr is scattered upward above the Ken- 
nicutt law for a given value of Nm, resulting in a broader 
deviation from the power-law relationship of the Kennicutt 
law near logA?Hi — 20. At higher values of Nm, Ssfr is en- 
hanced, but still within the range of the Kennicutt law. This 
diffe rence can be alleviated by a djusting the SF timescale 
fsee ISpringel fc Hernquist]|2003al ). however, we try to em- 
phasise the effect of metal cooing on SFR by keeping this 
parameter the same in all the runs. Further investigations 
of the Hi aspect of our simulations is beyond the scope of 
the present paper, and we will present the results elsewhere. 

Upon evaluatin g the effect of metal cooling, 
ISpringel fc HernauistJ l|2003bl ) considered two different 
gas phases where the cooling become important. One is the 
diffuse gas in galactic halos, which must radiate its thermal 
energy in order to collapse onto the high-density ISM, and 
the other is the multiphase star-forming gas. They argued 
in their § 5.2 that the cosmic SFR density at z > 6 would 
not be affected by the metal cooling very much, because 
at high-2 cooling is so efficient that the gas in diffuse 
atmospheres of halos cools nearly instantly, even without 
any metal cooling. They expected their model results to 
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Figure 10. Panel (a) highlights the redshift evolution of / gas in the N288L34 run from z = 5 to z = 1. Panel (b) shows the effect of 
numerical resolution on / ga s- The vertical dotted lines indicate the galaxy masses with 32 star particles for each run. The simulation 
results shift to lower masses as the resolution is increased. We plot the 'average' case in these figures. 
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Figure 11. Comparison of the N216L10 series with the empirical Kennicutt law. 



be largely independent of metal enrichment, because in 
this regime the evolution of SFR density is driven by 
the fast gravitational growth of the halo mass function. 
For the star-forming multiphase ISM, they argued that 
the parameters of their model (e.g. cold gas evaporation 
efficiency and gas consumption time-scale) can be adjusted 
such that the normalisation of the Kennicutt law can 
be maintained, yielding to first order unaltered model 
predictions. However, our simulations, in particular the 
comparison between the N216L10 and N216L10mc run, 
suggest that the metal cooling can enhance the SFR density 
even at z > 6 by ~ 20 — 30%. This is because the gas 
density becomes higher and temperature becomes lower 



(see Figure [5]:) by responding to the enhanced cooling rate 
by metals. 

In summary, we consider that our study serves to 
demonstrate the importance of metal enrichment and metal 
cooling in galaxy formation, beyond our original motivation 
to perform a simple numerical study. Our study highlights 
the role of metals at low-z and high-z to enhance the star 
formation in the Universe, and demonstrates that the metals 
increase the SF efficiency in the star-forming ISM, as well as 
enabling more IGM to accrete onto galaxies and fuel star for- 
mation. We expect that metal cooling would also alter some 
of the galaxy properties such as th e mass-metallicity rela- 
tionship (e.g.. IrTemonti et al.ll20o"3 ) and specific SFR (e.g., 
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iNoeske et ai1l2007f ). We plan to investigate these issues in 
the future, as well as alternative models for star formation 
and feedback. 
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